PROGRAM main
  USE global
  USE qVSCFintegrals
  USE forceconstants
  IMPLICIT NONE
  INTEGER:: mode
  REAL*8::lambda,nu
  Ndof=3
  nmr=3
  mode=1
  lambda=-0.1116d0
  nu=0.08414d0
  ALLOCATE(maxbasis(Ndof),hrmfreq(Ndof))
  maxbasis=4
  ALLOCATE(coef(nDOF,maxbasis(mode)))

  coef=0.0
  coef(:,1)=1.d0
  print*,coef
  CALL init_qff()
  CALL read_qff()
  CALL make_quartic()
  !CALL make_harmonic()
  !   CALL make_zero()
  !  Hii(1)=0.5*(0.29375d0)
  !  Hii(2)=0.5*(2.12581d0)
  !  Ciij(2,1)=lambda
  ! Ciii(1)=lambda*nu
  CALL make_units()
  ! CALL write_qff()
  hrmfreq(1)=1625.9276d0
  hrmfreq(2)=3799.6666d0
  hrmfreq(3)=3923.3194d0
  PRINT*,hrmfreq
  !  hrmfreq(1)=dsqrt(0.29375d0)
  !  hrmfreq(2)=dsqrt(2.12581d0)
  CALL vscf()
 ! PRINT*,'uzero',uzero(1),0.5*SQRT(uzero(1))
  !PRINT*,'utwo',utwo(1),SQRT(utwo(1))
  PRINT*,'nMR=',nMR
  PRINT*,'nDOF',nDOF
  PRINT*,'HO <000>=',SUM(hrmfreq)*0.5d0
  STOP
END PROGRAM main

